Assignment 3

Due date: Thursday 10/8, 6pm

This assignment will contain two parts:

  1. Exploring evictions and code violations in Philadelphia
  2. Comparing the NDVI in Philadelphia

Part 1: Exploring Evictions and Code Violations in Philadelphia

In this assignment, we'll explore spatial trends evictions in Philadelphia using data from the Eviction Lab and building code violations using data from OpenDataPhilly.

We'll be exploring the idea that evictions can occur as retaliation against renters for reporting code violations. Spatial correlations between evictions and code violations from the City's Licenses and Inspections department can offer some insight into this question.

A couple of interesting background readings:

1.1 Explore Eviction Lab Data

The Eviction Lab built the first national database for evictions. If you aren't familiar with the project, you can explore their website: https://evictionlab.org/

1.1.1 Read data using geopandas

The first step is to read the eviction data by census tract using geopandas. The data for all of Pennsylvania by census tract can be downloaded in a GeoJSON format using the following url:

https://eviction-lab-data-downloads.s3.amazonaws.com/PA/tracts.geojson

A browser-friendly version of the data is available here: https://data-downloads.evictionlab.org/

In [1]:
import numpy as np
from matplotlib import pyplot as plt
import pandas as pd
import geopandas as gpd

import hvplot.pandas
import holoviews as hv
hv.extension("bokeh")
In [2]:
#read eviction file
eviction = gpd.read_file('F:/MUSA/MUSA550/assignment-3/assignment-3-master/data/tracts.geojson')
eviction.head()
Out[2]:
GEOID west south east north n pl p-00 pr-00 roh-00 ... pm-16 po-16 ef-16 e-16 er-16 efr-16 lf-16 imputed-16 subbed-16 geometry
0 42003412002 -80.1243 40.5422 -80.0640 40.5890 4120.02 Allegheny County, Pennsylvania 4748.59 0.88 58.0 ... 0.00 0.0 0.0 0.0 0.00 0.00 1.0 0.0 1.0 MULTIPOLYGON (((-80.06670 40.58401, -80.06655 ...
1 42003413100 -80.0681 40.5850 -79.9906 40.6143 4131 Allegheny County, Pennsylvania 6771.01 3.47 729.0 ... 1.59 0.0 12.0 2.0 0.27 1.62 1.0 0.0 1.0 MULTIPOLYGON (((-80.06806 40.61254, -80.05452 ...
2 42003413300 -80.0657 40.5527 -80.0210 40.5721 4133 Allegheny County, Pennsylvania 5044.59 2.99 119.0 ... 0.95 0.0 4.0 1.0 0.49 1.96 1.0 0.0 1.0 MULTIPOLYGON (((-80.03822 40.55349, -80.04368 ...
3 42003416000 -79.8113 40.5440 -79.7637 40.5630 4160 Allegheny County, Pennsylvania 1775.93 4.99 121.0 ... 0.55 0.0 1.0 1.0 0.65 0.65 1.0 0.0 1.0 MULTIPOLYGON (((-79.76595 40.55092, -79.76542 ...
4 42003417200 -79.7948 40.5341 -79.7642 40.5443 4172 Allegheny County, Pennsylvania 1428.03 11.95 321.0 ... 0.00 0.0 7.0 3.0 0.82 1.90 1.0 0.0 1.0 MULTIPOLYGON (((-79.77114 40.54415, -79.76417 ...

5 rows × 399 columns

1.1.2 Explore and trim the data

We will need to trim data to Philadelphia only. Take a look at the data dictionary for the descriptions of the various columns: https://eviction-lab-data-downloads.s3.amazonaws.com/DATA_DICTIONARY.txt

Note: the column names are shortened — see the end of the above file for the abbreviations. The numbers at the end of the columns indicate the years. For example, e-16 is the number of evictions in 2016.

Take a look at the individual columns and trim to census tracts in Philadelphia. (Hint: Philadelphia is both a city and a county).

In [3]:
#trim to census tracts in Philadelphia
evict_philly = eviction.loc[eviction["pl"] == "Philadelphia County, Pennsylvania"].copy()
evict_philly.head()
Out[3]:
GEOID west south east north n pl p-00 pr-00 roh-00 ... pm-16 po-16 ef-16 e-16 er-16 efr-16 lf-16 imputed-16 subbed-16 geometry
435 42101000100 -75.1523 39.9481 -75.1415 39.9569 1 Philadelphia County, Pennsylvania 2646.71 9.26 1347.0 ... 2.49 0.00 25.0 16.0 0.93 1.45 0.0 0.0 1.0 MULTIPOLYGON (((-75.14161 39.95549, -75.14163 ...
436 42101000200 -75.1631 39.9529 -75.1511 39.9578 2 Philadelphia County, Pennsylvania 1362.00 56.42 374.0 ... 2.27 0.00 11.0 8.0 0.95 1.30 0.0 0.0 1.0 MULTIPOLYGON (((-75.15122 39.95686, -75.15167 ...
437 42101000300 -75.1798 39.9544 -75.1623 39.9599 3 Philadelphia County, Pennsylvania 2570.00 12.16 861.0 ... 1.76 0.00 26.0 14.0 0.73 1.35 0.0 0.0 1.0 MULTIPOLYGON (((-75.16234 39.95782, -75.16237 ...
438 42101000801 -75.1833 39.9486 -75.1773 39.9515 8.01 Philadelphia County, Pennsylvania 1478.00 14.40 810.0 ... 1.42 3.78 13.0 4.0 0.51 1.64 0.0 0.0 1.0 MULTIPOLYGON (((-75.17732 39.95096, -75.17784 ...
439 42101000804 -75.1712 39.9470 -75.1643 39.9501 8.04 Philadelphia County, Pennsylvania 3301.00 14.40 2058.0 ... 0.19 0.35 22.0 7.0 0.33 1.04 0.0 0.0 1.0 MULTIPOLYGON (((-75.17118 39.94778, -75.17102 ...

5 rows × 399 columns

1.1.3 Transform from wide to tidy format

For this assignment, we are interested in the number of evictions by census tract for various years. Right now, each year has it's own column, so it will be easiest to transform to a tidy format.

Use the pd.melt() function to transform the eviction data into tidy format, using the number of evictions from 2003 to 2016.

The tidy data frame should have four columns: GEOID, geometry, a column holding the number of evictions, and a column telling you what the name of the original column was for that value.

Hints:

  • You'll want to specify the GEOID and geometry columns as the id_vars. This will keep track of the census tract information.
  • You should specify the names of the columns holding the number of evictions as the value_vars.
  • You can generate a list of this column names using [Python's string formatting]:(https://docs.python.org/3.7/library/string.html#format-examples)
    value_vars = ['e-{:02d}'.format(x) for x in range(3, 17)]
    
In [4]:
#trim neccessary columns
value_vars = ['e-{:02d}'.format(x) for x in range(3, 17)]
value_vars.extend(["GEOID","geometry"])
evict_t = evict_philly[value_vars]

#transform the eviction data into tidy format
evict_p = pd.melt(evict_t, id_vars=["GEOID","geometry"], value_name="Number of evictions",var_name='Year')

#rename the column
key = ['e-{:02d}'.format(x) for x in range(3, 17)]
year = [i for i in range(2003,2017)]
replace_year = dict(zip(key,year))
evict_p['Year'].replace(replace_year,inplace=True)

evict_p.head()
Out[4]:
GEOID geometry Year Number of evictions
0 42101000100 MULTIPOLYGON (((-75.14161 39.95549, -75.14163 ... 2003 21.0
1 42101000200 MULTIPOLYGON (((-75.15122 39.95686, -75.15167 ... 2003 3.0
2 42101000300 MULTIPOLYGON (((-75.16234 39.95782, -75.16237 ... 2003 17.0
3 42101000801 MULTIPOLYGON (((-75.17732 39.95096, -75.17784 ... 2003 13.0
4 42101000804 MULTIPOLYGON (((-75.17118 39.94778, -75.17102 ... 2003 21.0

1.1.4 Plot the total number of evictions per year from 2003 to 2016

Use hvplot to plot the total number of evictions from 2003 to 2016. You will first need to perform a group by operation and sum up the total number of evictions for all census tracts, and then use hvplot() to make your plot.

You can use any type of hvplot chart you'd like to show the trend in number of evictions over time.

In [5]:
#group by year and count number of evictions
evict_sum = evict_p.groupby(['Year'])['Number of evictions'].sum() 
In [6]:
plot1 = evict_sum.hvplot(kind='line')
plot1
Out[6]:

1.1.5 The number of evictions across Philadelphia

Our tidy data frame is still a GeoDataFrame with a geometry column, so we can visualize the number of evictions for all census tracts.

Use hvplot() to generate a choropleth showing the number of evictions for a specified year, with a widget dropdown to select a given year (or variable name, e.g., e-16, e-15, etc).

Hints

  • You'll need to use the groupby keyword to tell hvplot to make a series of maps, with a widget to select between them.
  • You will need to specify dynamic=False as a keyword argument to the hvplot() function.
  • Be sure to specify a width and height that makes your output map (roughly) square to limit distortions
In [7]:
#plot number of evictions across Philadelphia 2003-2016
evict_p.hvplot(c='Number of evictions',
                     geo=True,
                     frame_width=500, 
                     frame_height=500,
                     groupby="Year",
                     dynamic=False,
                     cmap='viridis',
                     title='The number of evictions across Philadelphia')
F:\anaconda\envs\musa-550-fall-2020\lib\site-packages\holoviews\plotting\util.py:685: MatplotlibDeprecationWarning: The global colormaps dictionary is no longer considered public API.
  [cmap for cmap in cm.cmap_d if not
F:\anaconda\envs\musa-550-fall-2020\lib\site-packages\holoviews\plotting\util.py:685: MatplotlibDeprecationWarning: The global colormaps dictionary is no longer considered public API.
  [cmap for cmap in cm.cmap_d if not
F:\anaconda\envs\musa-550-fall-2020\lib\site-packages\holoviews\plotting\util.py:685: MatplotlibDeprecationWarning: The global colormaps dictionary is no longer considered public API.
  [cmap for cmap in cm.cmap_d if not
F:\anaconda\envs\musa-550-fall-2020\lib\site-packages\holoviews\plotting\util.py:685: MatplotlibDeprecationWarning: The global colormaps dictionary is no longer considered public API.
  [cmap for cmap in cm.cmap_d if not
F:\anaconda\envs\musa-550-fall-2020\lib\site-packages\holoviews\plotting\util.py:685: MatplotlibDeprecationWarning: The global colormaps dictionary is no longer considered public API.
  [cmap for cmap in cm.cmap_d if not
F:\anaconda\envs\musa-550-fall-2020\lib\site-packages\holoviews\plotting\util.py:685: MatplotlibDeprecationWarning: The global colormaps dictionary is no longer considered public API.
  [cmap for cmap in cm.cmap_d if not
F:\anaconda\envs\musa-550-fall-2020\lib\site-packages\holoviews\plotting\util.py:685: MatplotlibDeprecationWarning: The global colormaps dictionary is no longer considered public API.
  [cmap for cmap in cm.cmap_d if not
F:\anaconda\envs\musa-550-fall-2020\lib\site-packages\holoviews\plotting\util.py:685: MatplotlibDeprecationWarning: The global colormaps dictionary is no longer considered public API.
  [cmap for cmap in cm.cmap_d if not
F:\anaconda\envs\musa-550-fall-2020\lib\site-packages\holoviews\plotting\util.py:685: MatplotlibDeprecationWarning: The global colormaps dictionary is no longer considered public API.
  [cmap for cmap in cm.cmap_d if not
F:\anaconda\envs\musa-550-fall-2020\lib\site-packages\holoviews\plotting\util.py:685: MatplotlibDeprecationWarning: The global colormaps dictionary is no longer considered public API.
  [cmap for cmap in cm.cmap_d if not
F:\anaconda\envs\musa-550-fall-2020\lib\site-packages\holoviews\plotting\util.py:685: MatplotlibDeprecationWarning: The global colormaps dictionary is no longer considered public API.
  [cmap for cmap in cm.cmap_d if not
F:\anaconda\envs\musa-550-fall-2020\lib\site-packages\holoviews\plotting\util.py:685: MatplotlibDeprecationWarning: The global colormaps dictionary is no longer considered public API.
  [cmap for cmap in cm.cmap_d if not
F:\anaconda\envs\musa-550-fall-2020\lib\site-packages\holoviews\plotting\util.py:685: MatplotlibDeprecationWarning: The global colormaps dictionary is no longer considered public API.
  [cmap for cmap in cm.cmap_d if not
F:\anaconda\envs\musa-550-fall-2020\lib\site-packages\holoviews\plotting\util.py:685: MatplotlibDeprecationWarning: The global colormaps dictionary is no longer considered public API.
  [cmap for cmap in cm.cmap_d if not
F:\anaconda\envs\musa-550-fall-2020\lib\site-packages\holoviews\plotting\util.py:685: MatplotlibDeprecationWarning: The global colormaps dictionary is no longer considered public API.
  [cmap for cmap in cm.cmap_d if not
F:\anaconda\envs\musa-550-fall-2020\lib\site-packages\holoviews\plotting\util.py:685: MatplotlibDeprecationWarning: The global colormaps dictionary is no longer considered public API.
  [cmap for cmap in cm.cmap_d if not
F:\anaconda\envs\musa-550-fall-2020\lib\site-packages\holoviews\plotting\util.py:685: MatplotlibDeprecationWarning: The global colormaps dictionary is no longer considered public API.
  [cmap for cmap in cm.cmap_d if not
F:\anaconda\envs\musa-550-fall-2020\lib\site-packages\holoviews\plotting\util.py:685: MatplotlibDeprecationWarning: The global colormaps dictionary is no longer considered public API.
  [cmap for cmap in cm.cmap_d if not
Out[7]:

1.2 Code Violations in Philadelphia

Next, we'll explore data for code violations from the Licenses and Inspections Department of Philadelphia to look for potential correlations with the number of evictions.

1.2.1 Load data from 2012 to 2016

L+I violation data for years including 2012 through 2016 (inclusive) is provided in a CSV format in the "data/" folder.

Load the data using pandas and convert to a GeoDataFrame.

In [8]:
#read csv data and add geometry
violation = pd.read_csv('F:/MUSA/MUSA550/assignment-3/assignment-3-master/data/li_violations.csv')
violation = violation.dropna(subset=['lat', 'lng']) 
violation['Coordinates'] = gpd.points_from_xy(violation['lng'], violation['lat'])
violation = gpd.GeoDataFrame(violation, 
                            geometry="Coordinates", 
                            crs="EPSG:4326")
violation.head()
Out[8]:
lat lng violationdescription Coordinates
0 40.050526 -75.126076 CLIP VIOLATION NOTICE POINT (-75.12608 40.05053)
1 40.050593 -75.126578 LICENSE-CHANGE OF ADDRESS POINT (-75.12658 40.05059)
2 40.050593 -75.126578 LICENSE-RES SFD/2FD POINT (-75.12658 40.05059)
3 39.991994 -75.128895 EXT A-CLEAN WEEDS/PLANTS POINT (-75.12889 39.99199)
4 40.023260 -75.164848 EXT A-VACANT LOT CLEAN/MAINTAI POINT (-75.16485 40.02326)

1.2.2 Trim to specific violation types

There are many different types of code violations (running the nunique() function on the violationdescription column will extract all of the unique ones). More information on different types of violations can be found on the City's website.

Below, I've selected 15 types of violations that deal with property maintenance and licensing issues. We'll focus on these violations. The goal is to see if these kinds of violations are correlated spatially with the number of evictions in a given area.

Use the list of violations given to trim your data set to only include these types.

In [9]:
violation_types = [
    "INT-PLMBG MAINT FIXTURES-RES",
    "INT S-CEILING REPAIR/MAINT SAN",
    "PLUMBING SYSTEMS-GENERAL",
    "CO DETECTOR NEEDED",
    "INTERIOR SURFACES",
    "EXT S-ROOF REPAIR",
    "ELEC-RECEPTABLE DEFECTIVE-RES",
    "INT S-FLOOR REPAIR",
    "DRAINAGE-MAIN DRAIN REPAIR-RES",
    "DRAINAGE-DOWNSPOUT REPR/REPLC",
    "LIGHT FIXTURE DEFECTIVE-RES",
    "LICENSE-RES SFD/2FD",
    "ELECTRICAL -HAZARD",
    "VACANT PROPERTIES-GENERAL",
    "INT-PLMBG FIXTURES-RES",
]
In [10]:
# Trim to specific violation types
violation_t = violation.loc[violation["violationdescription"].isin(violation_types)].copy()
violation_t.head()
Out[10]:
lat lng violationdescription Coordinates
2 40.050593 -75.126578 LICENSE-RES SFD/2FD POINT (-75.12658 40.05059)
25 40.022406 -75.121872 EXT S-ROOF REPAIR POINT (-75.12187 40.02241)
30 40.023237 -75.121726 CO DETECTOR NEEDED POINT (-75.12173 40.02324)
31 40.023397 -75.122241 INT S-CEILING REPAIR/MAINT SAN POINT (-75.12224 40.02340)
34 40.023773 -75.121603 INT S-FLOOR REPAIR POINT (-75.12160 40.02377)

1.2.3 Make a hex bin map

The code violation data is point data. We can get a quick look at the geographic distribution using matplotlib and the hexbin() function. Make a hex bin map of the code violations and overlay the census tract outlines.

Hints:

  • The eviction data from part 1 was by census tract, so the census tract geometries are available as part of that GeoDataFrame. You can use it to overlay the census tracts on your hex bin map.
  • Make sure you convert your GeoDataFrame to a CRS that's better for visualization than plain old 4326.
In [11]:
# convert to same CRS
evict_3857 = evict_t[['geometry','GEOID']].to_crs(epsg=3857)
violation_3857 = violation_t.to_crs(epsg=3857)

fig, ax = plt.subplots(figsize=(18, 12))


# Extract out the x/y coordindates of the Point objects
xcoords = violation_3857.Coordinates.x
ycoords = violation_3857.Coordinates.y

# Plot a hexbin chart
hex_vals = ax.hexbin(xcoords, ycoords, gridsize=50)

# Add the tracts
evict_3857.plot(ax=ax, facecolor="none", edgecolor="white",alpha=0.75, linewidth=0.25)

# Get the limits of the GeoDataFrame
xmin, ymin, xmax, ymax = evict_3857.total_bounds

# Set the xlims and ylims
ax.set_xlim(xmin, xmax)
ax.set_ylim(ymin, ymax)

# add a colorbar and format
fig.colorbar(hex_vals, ax=ax)
ax.set_axis_off()
ax.set_aspect("equal")
ax.set_title("Violation in Philadelphia 2012-2016")
Out[11]:
Text(0.5, 1.0, 'Violation in Philadelphia 2012-2016')

1.2.4 Spatially join data sets

To do a census tract comparison to our eviction data, we need to find which census tract each of the code violations falls into. Use the geopandas.sjoin() function to do just that.

Hints

  • You can re-use your eviction data frame, but you will only need the geometry column (specifying census tract polygons) and the GEOID column (specifying the name of each census tract).
  • Make sure both data frames have the same CRS before joining them together!
In [12]:
#Spatially join data sets
joined = gpd.sjoin(violation_3857, evict_3857, op='within', how='right')
joined['GEOID'] = joined['GEOID'].astype('category')
len(joined['GEOID'].unique())
Out[12]:
384

1.2.5 Calculate the number of violations by type per census tract

Next, we'll want to find the number of violations (for each kind) per census tract. You should group the data frame by violation type and census tract name.

The result of this step should be a data frame with three columns: violationdescription, GEOID, and N, where N is the number of violations of that kind in the specified census tract.

Optional: to make prettier plots

Some census tracts won't have any violations, and they won't be included when we do the above calculation. However, there is a trick to set the values for those census tracts to be zero. After you calculate the sizes of each violation/census tract group, you can run:

N = N.unstack(fill_value=0).stack().reset_index(name='N')

where N gives the total size of each of the groups, specified by violation type and census tract name.

See this StackOverflow post for more details.

This part is optional, but will make the resulting maps a bit prettier.

In [13]:
vio_g = joined.groupby(['GEOID','violationdescription']).size().unstack(fill_value=0).stack().reset_index(name=0).rename(columns={0:'N'})
len(vio_g['GEOID'].unique())
Out[13]:
384

1.2.6 Merge with census tracts geometries

We now have the number of violations of different types per census tract specified as a regular DataFrame. You can now merge it with the census tract geometries (from your eviction data GeoDataFrame) to create a GeoDataFrame.

Hints

  • Use pandas.merge() and specify the on keyword to be the column holding census tract names.
  • Make sure the result of the merge operation is a GeoDataFrame — you will want the GeoDataFrame holding census tract geometries to be the first argument of the pandas.merge() function.
In [14]:
total = evict_3857.merge(vio_g, on='GEOID')
total.head()
Out[14]:
geometry GEOID violationdescription N
0 MULTIPOLYGON (((-8364725.429 4859476.459, -836... 42101000100 CO DETECTOR NEEDED 0
1 MULTIPOLYGON (((-8364725.429 4859476.459, -836... 42101000100 DRAINAGE-DOWNSPOUT REPR/REPLC 6
2 MULTIPOLYGON (((-8364725.429 4859476.459, -836... 42101000100 DRAINAGE-MAIN DRAIN REPAIR-RES 0
3 MULTIPOLYGON (((-8364725.429 4859476.459, -836... 42101000100 ELEC-RECEPTABLE DEFECTIVE-RES 0
4 MULTIPOLYGON (((-8364725.429 4859476.459, -836... 42101000100 ELECTRICAL -HAZARD 1

1.2.7 Interactive choropleths for each violation type

Now, we can use hvplot() to create an interactive choropleth for each violation type and add a widget to specify different violation types.

Hints

  • You'll need to use the groupby keyword to tell hvplot to make a series of maps, with a widget to select different violation types.
  • You will need to specify dynamic=False as a keyword argument to the hvplot() function.
  • Be sure to specify a width and height that makes your output map (roughly) square to limit distortions
In [15]:
#plot interactive choropleths for each violation type
total=total.to_crs(epsg=4326)

total.hvplot(c='N',
            geo=True,
            frame_width=500, 
            frame_height=500,
            groupby="violationdescription",
            dynamic=False,
            cmap='viridis')
F:\anaconda\envs\musa-550-fall-2020\lib\site-packages\holoviews\plotting\util.py:685: MatplotlibDeprecationWarning: The global colormaps dictionary is no longer considered public API.
  [cmap for cmap in cm.cmap_d if not
F:\anaconda\envs\musa-550-fall-2020\lib\site-packages\holoviews\plotting\util.py:685: MatplotlibDeprecationWarning: The global colormaps dictionary is no longer considered public API.
  [cmap for cmap in cm.cmap_d if not
F:\anaconda\envs\musa-550-fall-2020\lib\site-packages\holoviews\plotting\util.py:685: MatplotlibDeprecationWarning: The global colormaps dictionary is no longer considered public API.
  [cmap for cmap in cm.cmap_d if not
F:\anaconda\envs\musa-550-fall-2020\lib\site-packages\holoviews\plotting\util.py:685: MatplotlibDeprecationWarning: The global colormaps dictionary is no longer considered public API.
  [cmap for cmap in cm.cmap_d if not
F:\anaconda\envs\musa-550-fall-2020\lib\site-packages\holoviews\plotting\util.py:685: MatplotlibDeprecationWarning: The global colormaps dictionary is no longer considered public API.
  [cmap for cmap in cm.cmap_d if not
F:\anaconda\envs\musa-550-fall-2020\lib\site-packages\holoviews\plotting\util.py:685: MatplotlibDeprecationWarning: The global colormaps dictionary is no longer considered public API.
  [cmap for cmap in cm.cmap_d if not
F:\anaconda\envs\musa-550-fall-2020\lib\site-packages\holoviews\plotting\util.py:685: MatplotlibDeprecationWarning: The global colormaps dictionary is no longer considered public API.
  [cmap for cmap in cm.cmap_d if not
F:\anaconda\envs\musa-550-fall-2020\lib\site-packages\holoviews\plotting\util.py:685: MatplotlibDeprecationWarning: The global colormaps dictionary is no longer considered public API.
  [cmap for cmap in cm.cmap_d if not
F:\anaconda\envs\musa-550-fall-2020\lib\site-packages\holoviews\plotting\util.py:685: MatplotlibDeprecationWarning: The global colormaps dictionary is no longer considered public API.
  [cmap for cmap in cm.cmap_d if not
F:\anaconda\envs\musa-550-fall-2020\lib\site-packages\holoviews\plotting\util.py:685: MatplotlibDeprecationWarning: The global colormaps dictionary is no longer considered public API.
  [cmap for cmap in cm.cmap_d if not
F:\anaconda\envs\musa-550-fall-2020\lib\site-packages\holoviews\plotting\util.py:685: MatplotlibDeprecationWarning: The global colormaps dictionary is no longer considered public API.
  [cmap for cmap in cm.cmap_d if not
F:\anaconda\envs\musa-550-fall-2020\lib\site-packages\holoviews\plotting\util.py:685: MatplotlibDeprecationWarning: The global colormaps dictionary is no longer considered public API.
  [cmap for cmap in cm.cmap_d if not
F:\anaconda\envs\musa-550-fall-2020\lib\site-packages\holoviews\plotting\util.py:685: MatplotlibDeprecationWarning: The global colormaps dictionary is no longer considered public API.
  [cmap for cmap in cm.cmap_d if not
F:\anaconda\envs\musa-550-fall-2020\lib\site-packages\holoviews\plotting\util.py:685: MatplotlibDeprecationWarning: The global colormaps dictionary is no longer considered public API.
  [cmap for cmap in cm.cmap_d if not
F:\anaconda\envs\musa-550-fall-2020\lib\site-packages\holoviews\plotting\util.py:685: MatplotlibDeprecationWarning: The global colormaps dictionary is no longer considered public API.
  [cmap for cmap in cm.cmap_d if not
F:\anaconda\envs\musa-550-fall-2020\lib\site-packages\holoviews\plotting\util.py:685: MatplotlibDeprecationWarning: The global colormaps dictionary is no longer considered public API.
  [cmap for cmap in cm.cmap_d if not
F:\anaconda\envs\musa-550-fall-2020\lib\site-packages\holoviews\plotting\util.py:685: MatplotlibDeprecationWarning: The global colormaps dictionary is no longer considered public API.
  [cmap for cmap in cm.cmap_d if not
F:\anaconda\envs\musa-550-fall-2020\lib\site-packages\holoviews\plotting\util.py:685: MatplotlibDeprecationWarning: The global colormaps dictionary is no longer considered public API.
  [cmap for cmap in cm.cmap_d if not
F:\anaconda\envs\musa-550-fall-2020\lib\site-packages\holoviews\plotting\util.py:685: MatplotlibDeprecationWarning: The global colormaps dictionary is no longer considered public API.
  [cmap for cmap in cm.cmap_d if not
Out[15]:

1.3. A side-by-side comparison

From the interactive maps of evictions and violations, you should notice a lot of spatial overlap.

As a final step, we'll make a side-by-side comparison to better show the spatial correlations. This will involve a few steps:

  1. Trim the data frame plotted in section 1.1.5 to only include evictions from 2016.
  2. Trim the data frame plotted in section 1.2.7 to only include a single violation type (pick whichever one you want!).
  3. Use hvplot() to make two interactive choropleth maps, one for the data from step 1. and one for the data in step 2.
  4. Show these two plots side by side (one row and 2 columns) using the syntax for combining charts.

Note: since we selected a single year and violation type, you won't need to use the groupby= keyword here.

In [16]:
#Trim the data frame plotted in section 1.1.5 to only include evictions from 2016.
evict_2016 = evict_p.loc[evict_p['Year']==2016].reset_index()
In [17]:
#Trim the data frame plotted in section 1.2.7 to only include a single violation type (pick whichever one you want!)
vio_vp = total.loc[total['violationdescription']=="VACANT PROPERTIES-GENERAL"].reset_index()
In [18]:
#Use hvplot() to make two interactive choropleth maps, one for the data from step 1. and one for the data in step 2.
plot1=evict_2016.hvplot(c='Number of evictions',
            geo=True,
            frame_width=500, 
            frame_height=500,
            dynamic=False,
            cmap='viridis')

plot1
F:\anaconda\envs\musa-550-fall-2020\lib\site-packages\holoviews\plotting\util.py:685: MatplotlibDeprecationWarning: The global colormaps dictionary is no longer considered public API.
  [cmap for cmap in cm.cmap_d if not
F:\anaconda\envs\musa-550-fall-2020\lib\site-packages\holoviews\plotting\util.py:685: MatplotlibDeprecationWarning: The global colormaps dictionary is no longer considered public API.
  [cmap for cmap in cm.cmap_d if not
Out[18]:
In [19]:
plot2=vio_vp.hvplot(c='N',
            geo=True,
            frame_width=500, 
            frame_height=500,
            dynamic=False,
            cmap='viridis')

plot2
F:\anaconda\envs\musa-550-fall-2020\lib\site-packages\holoviews\plotting\util.py:685: MatplotlibDeprecationWarning: The global colormaps dictionary is no longer considered public API.
  [cmap for cmap in cm.cmap_d if not
F:\anaconda\envs\musa-550-fall-2020\lib\site-packages\holoviews\plotting\util.py:685: MatplotlibDeprecationWarning: The global colormaps dictionary is no longer considered public API.
  [cmap for cmap in cm.cmap_d if not
Out[19]:
In [20]:
#Show these two plots side by side (one row and 2 columns) using the syntax for combining charts.
plot1=plot1.relabel('Number of evictions 2016')
plot2=plot2.relabel('Vacant properties general')

combined = plot1 + plot2 
In [21]:
combined.cols(2)
F:\anaconda\envs\musa-550-fall-2020\lib\site-packages\holoviews\plotting\util.py:685: MatplotlibDeprecationWarning: The global colormaps dictionary is no longer considered public API.
  [cmap for cmap in cm.cmap_d if not
F:\anaconda\envs\musa-550-fall-2020\lib\site-packages\holoviews\plotting\util.py:685: MatplotlibDeprecationWarning: The global colormaps dictionary is no longer considered public API.
  [cmap for cmap in cm.cmap_d if not
Out[21]:

1.4. Extra Credit

Identify the 20 most common types of violations within the time period of 2012 to 2016 and create a set of interactive choropleths similar to what was done in section 1.2.7.

Use this set of maps to identify 3 types of violations that don't seem to have much spatial overlap with the number of evictions in the City.

In [22]:
#trim the table with top20 violations
top20 = violation.groupby('violationdescription').size()
top20 = top20.sort_values(ascending=False)
top20 = top20.iloc[:20].reset_index()
list20 = top20['violationdescription'].values.tolist()
vio20 = violation.loc[violation["violationdescription"].isin(list20)].copy()
list20
Out[22]:
['CLIP VIOLATION NOTICE',
 'EXT A-VACANT LOT CLEAN/MAINTAI',
 'HIGH WEEDS-CUT',
 'LICENSE-VAC RES BLDG',
 'VACANT PROP STANDARD',
 'RUBBISH/GARBAGE EXTERIOR-OWNER',
 'EXT A-CLEAN RUBBISH/GARBAGE',
 'LICENSE-RES SFD/2FD',
 'EXT A-CLEAN WEEDS/PLANTS',
 'LICENSE-RES GENERAL',
 'VACANT BLDG UNSECURED COUNT',
 'INT S-CEILING REPAIR/MAINT SAN',
 'VIOL C&I MESSAGE',
 'CO DETECTOR NEEDED',
 'ANNUAL CERT FIRE ALARM',
 'LICENSE - RENTAL PROPERTY',
 'VAC PROP REPLAC WIN/DRS 80%',
 'SD-REQD EXIST GROUP R',
 'PERM Z- NEW USE',
 'INT S-WALLS REPAIR/MAINT SANI']
In [23]:
#spatial join
vio20_3857=vio20.to_crs(epsg=3857)
joined1 = gpd.sjoin(vio20_3857, evict_3857, op='within', how='right')
joined1['GEOID'] = joined1['GEOID'].astype('category')
len(joined1['GEOID'].unique())
Out[23]:
384
In [24]:
#Calculate the number of violations by type per census tract¶
vio20_g = joined1.groupby(['GEOID','violationdescription']).size().unstack(fill_value=0).stack().reset_index(name=0).rename(columns={0:'N'})
len(vio20_g['GEOID'].unique())
Out[24]:
384
In [25]:
#merge the violation data with tracts
total1 = evict_3857.merge(vio20_g, on='GEOID')
total.head()
Out[25]:
geometry GEOID violationdescription N
0 MULTIPOLYGON (((-75.14161 39.95549, -75.14163 ... 42101000100 CO DETECTOR NEEDED 0
1 MULTIPOLYGON (((-75.14161 39.95549, -75.14163 ... 42101000100 DRAINAGE-DOWNSPOUT REPR/REPLC 6
2 MULTIPOLYGON (((-75.14161 39.95549, -75.14163 ... 42101000100 DRAINAGE-MAIN DRAIN REPAIR-RES 0
3 MULTIPOLYGON (((-75.14161 39.95549, -75.14163 ... 42101000100 ELEC-RECEPTABLE DEFECTIVE-RES 0
4 MULTIPOLYGON (((-75.14161 39.95549, -75.14163 ... 42101000100 ELECTRICAL -HAZARD 1
In [26]:
total1=total1.to_crs(epsg=4326)

total1.hvplot(c='N',
            geo=True,
            frame_width=500, 
            frame_height=500,
            groupby="violationdescription",
            dynamic=False,
            cmap='viridis')
F:\anaconda\envs\musa-550-fall-2020\lib\site-packages\holoviews\plotting\util.py:685: MatplotlibDeprecationWarning: The global colormaps dictionary is no longer considered public API.
  [cmap for cmap in cm.cmap_d if not
F:\anaconda\envs\musa-550-fall-2020\lib\site-packages\holoviews\plotting\util.py:685: MatplotlibDeprecationWarning: The global colormaps dictionary is no longer considered public API.
  [cmap for cmap in cm.cmap_d if not
F:\anaconda\envs\musa-550-fall-2020\lib\site-packages\holoviews\plotting\util.py:685: MatplotlibDeprecationWarning: The global colormaps dictionary is no longer considered public API.
  [cmap for cmap in cm.cmap_d if not
F:\anaconda\envs\musa-550-fall-2020\lib\site-packages\holoviews\plotting\util.py:685: MatplotlibDeprecationWarning: The global colormaps dictionary is no longer considered public API.
  [cmap for cmap in cm.cmap_d if not
F:\anaconda\envs\musa-550-fall-2020\lib\site-packages\holoviews\plotting\util.py:685: MatplotlibDeprecationWarning: The global colormaps dictionary is no longer considered public API.
  [cmap for cmap in cm.cmap_d if not
F:\anaconda\envs\musa-550-fall-2020\lib\site-packages\holoviews\plotting\util.py:685: MatplotlibDeprecationWarning: The global colormaps dictionary is no longer considered public API.
  [cmap for cmap in cm.cmap_d if not
F:\anaconda\envs\musa-550-fall-2020\lib\site-packages\holoviews\plotting\util.py:685: MatplotlibDeprecationWarning: The global colormaps dictionary is no longer considered public API.
  [cmap for cmap in cm.cmap_d if not
F:\anaconda\envs\musa-550-fall-2020\lib\site-packages\holoviews\plotting\util.py:685: MatplotlibDeprecationWarning: The global colormaps dictionary is no longer considered public API.
  [cmap for cmap in cm.cmap_d if not
F:\anaconda\envs\musa-550-fall-2020\lib\site-packages\holoviews\plotting\util.py:685: MatplotlibDeprecationWarning: The global colormaps dictionary is no longer considered public API.
  [cmap for cmap in cm.cmap_d if not
F:\anaconda\envs\musa-550-fall-2020\lib\site-packages\holoviews\plotting\util.py:685: MatplotlibDeprecationWarning: The global colormaps dictionary is no longer considered public API.
  [cmap for cmap in cm.cmap_d if not
F:\anaconda\envs\musa-550-fall-2020\lib\site-packages\holoviews\plotting\util.py:685: MatplotlibDeprecationWarning: The global colormaps dictionary is no longer considered public API.
  [cmap for cmap in cm.cmap_d if not
F:\anaconda\envs\musa-550-fall-2020\lib\site-packages\holoviews\plotting\util.py:685: MatplotlibDeprecationWarning: The global colormaps dictionary is no longer considered public API.
  [cmap for cmap in cm.cmap_d if not
F:\anaconda\envs\musa-550-fall-2020\lib\site-packages\holoviews\plotting\util.py:685: MatplotlibDeprecationWarning: The global colormaps dictionary is no longer considered public API.
  [cmap for cmap in cm.cmap_d if not
F:\anaconda\envs\musa-550-fall-2020\lib\site-packages\holoviews\plotting\util.py:685: MatplotlibDeprecationWarning: The global colormaps dictionary is no longer considered public API.
  [cmap for cmap in cm.cmap_d if not
F:\anaconda\envs\musa-550-fall-2020\lib\site-packages\holoviews\plotting\util.py:685: MatplotlibDeprecationWarning: The global colormaps dictionary is no longer considered public API.
  [cmap for cmap in cm.cmap_d if not
F:\anaconda\envs\musa-550-fall-2020\lib\site-packages\holoviews\plotting\util.py:685: MatplotlibDeprecationWarning: The global colormaps dictionary is no longer considered public API.
  [cmap for cmap in cm.cmap_d if not
F:\anaconda\envs\musa-550-fall-2020\lib\site-packages\holoviews\plotting\util.py:685: MatplotlibDeprecationWarning: The global colormaps dictionary is no longer considered public API.
  [cmap for cmap in cm.cmap_d if not
F:\anaconda\envs\musa-550-fall-2020\lib\site-packages\holoviews\plotting\util.py:685: MatplotlibDeprecationWarning: The global colormaps dictionary is no longer considered public API.
  [cmap for cmap in cm.cmap_d if not
F:\anaconda\envs\musa-550-fall-2020\lib\site-packages\holoviews\plotting\util.py:685: MatplotlibDeprecationWarning: The global colormaps dictionary is no longer considered public API.
  [cmap for cmap in cm.cmap_d if not
F:\anaconda\envs\musa-550-fall-2020\lib\site-packages\holoviews\plotting\util.py:685: MatplotlibDeprecationWarning: The global colormaps dictionary is no longer considered public API.
  [cmap for cmap in cm.cmap_d if not
F:\anaconda\envs\musa-550-fall-2020\lib\site-packages\holoviews\plotting\util.py:685: MatplotlibDeprecationWarning: The global colormaps dictionary is no longer considered public API.
  [cmap for cmap in cm.cmap_d if not
F:\anaconda\envs\musa-550-fall-2020\lib\site-packages\holoviews\plotting\util.py:685: MatplotlibDeprecationWarning: The global colormaps dictionary is no longer considered public API.
  [cmap for cmap in cm.cmap_d if not
F:\anaconda\envs\musa-550-fall-2020\lib\site-packages\holoviews\plotting\util.py:685: MatplotlibDeprecationWarning: The global colormaps dictionary is no longer considered public API.
  [cmap for cmap in cm.cmap_d if not
F:\anaconda\envs\musa-550-fall-2020\lib\site-packages\holoviews\plotting\util.py:685: MatplotlibDeprecationWarning: The global colormaps dictionary is no longer considered public API.
  [cmap for cmap in cm.cmap_d if not
Out[26]:

3 types of violations that don't seem to have much spatial overlap with the number of evictions in the City

1.EXT A-CLEAN RUBBISH/GARBAGE

2.LICENSE-RES GENERAL

3.VIOL C&I MESSAGE

Part 2: Exploring the NDVI in Philadelphia

In this part, we'll explore the NDVI in Philadelphia a bit more. This part will include two parts:

  1. We'll compare the median NDVI within the city limits and the immediate suburbs
  2. We'll calculate the NDVI around street trees in the city.

2.1 Comparing the NDVI in the city and the suburbs

2.1.1 Load Landsat data for Philadelphia

Use rasterio to load the landsat data for Philadelphia (available in the "data/" folder)

In [27]:
import rasterio as rio
In [28]:
# Open the file
landsat = rio.open("F:/MUSA/MUSA550/assignment-3/assignment-3-master/data/landsat8_philly.tif")
# Read the file
data = landsat.read(1)

2.1.2 Separating the city from the suburbs

Create two polygon objects, one for the city limits and one for the suburbs. To calculate the suburbs polygon, we will take everything outside the city limits but still within the bounding box.

  • The city limits are available in the "data/" folder.
  • To calculate the suburbs polygon, the "envelope" attribute of the city limits geometry will be useful.
  • You can use geopandas' geometric manipulation functionality to calculate the suburbs polygon from the city limits polygon and the envelope polygon.
In [29]:
#city polygon
city_limits = gpd.read_file("F:/MUSA/MUSA550/assignment-3/assignment-3-master/data/City_Limits.geojson")
                            
#envelop polygon
box = city_limits.envelope

#suburb polygon
suburb = box.difference(city_limits)

2.1.3 Mask and calculate the NDVI for the city and the suburbs

Using the two polygons from the last section, use rasterio's mask functionality to create two masked arrays from the landsat data, one for the city and one for the suburbs.

For each masked array, calculate the NDVI.

In [30]:
from rasterio.mask import mask
import matplotlib.colors as mcolors
In [31]:
#convert CRS
city_limits = city_limits.to_crs(landsat.crs.data['init'])

#mask for city
city_masked, mask_transform = mask(
    dataset=landsat,
    shapes=city_limits.geometry,
    crop=True, 
    all_touched=True,  
    filled=False,  
)
In [32]:
#convert CRS
sub_limits = suburb.to_crs(landsat.crs.data['init'])

#mask for suburb
sub_masked, mask_transform = mask(
    dataset=landsat,
    shapes=sub_limits.geometry,
    crop=True, 
    all_touched=True,  
    filled=False,  
)
In [33]:
#NDVI founction
def calculate_NDVI(nir, red):
    """
    Calculate the NDVI from the NIR and red landsat bands
    """
    
    # Convert to floats
    nir = nir.astype(float)
    red = red.astype(float)
    
    # Get valid entries
    check = np.logical_and( red.mask == False, nir.mask == False )
    
    # Where the check is True, return the NDVI, else return NaN
    ndvi = np.where(check,  (nir - red ) / ( nir + red ), np.nan )
    return ndvi 
In [34]:
#calculate NDVI
NDVI_city = calculate_NDVI(city_masked[4], city_masked[3])
NDVI_sub = calculate_NDVI(sub_masked[4], sub_masked[3])
In [40]:
#plot NDVI in Philadelphia City
fig, ax = plt.subplots(figsize=(10,10))

# The extent of the data
landsat_extent = [
    landsat.bounds.left,
    landsat.bounds.right,
    landsat.bounds.bottom,
    landsat.bounds.top,
]

# Plot NDVI
img = ax.imshow(NDVI_city, extent=landsat_extent)

# Format and plot city limits
city_limits.plot(ax=ax, edgecolor='pink', facecolor='none', linewidth=4)
plt.colorbar(img)
ax.set_axis_off()
ax.set_title("NDVI in Philadelphia City", fontsize=18);
In [41]:
#plot NDVI in Philadelphia suburbs
fig, ax = plt.subplots(figsize=(10,10))

# Plot NDVI
img = ax.imshow(NDVI_sub, extent=landsat_extent)

# Format and plot city limits
sub_limits.plot(ax=ax, edgecolor='pink', facecolor='none', linewidth=4)
plt.colorbar(img)
ax.set_axis_off()
ax.set_title("NDVI in Philadelphia Suburbs", fontsize=18);

2.1.4 Calculate the median NDVI within the city and within the suburbs

  • Calculate the median value from your NDVI arrays for the city and suburbs
  • Numpy's nanmedian function will be useful for ignoring NaN elements
  • Print out the median values. Which has a higher NDVI: the city or suburbs?
In [42]:
MNDVI_city=np.nanmedian(NDVI_city)
MNDVI_sub=np.nanmedian(NDVI_sub)
print(MNDVI_city)
print(MNDVI_sub)
0.20268593532493442
0.37466958688920776

Suburbs have higher median NDVI

2.2 Calculating the NDVI for Philadelphia's street trees

2.2.1 Load the street tree data

The data is available in the "data/" folder. It has been downloaded from OpenDataPhilly. It contains the locations of abot 2,500 street trees in Philadelphia.

In [43]:
#read the file
tree = gpd.read_file('F:/MUSA/MUSA550/assignment-3/assignment-3-master/data/ppr_tree_canopy_points_2015.geojson')
tree.head()
#convert CRS
tree = tree.to_crs(landsat.crs.data['init'])

2.2.2 Calculate the NDVI values at the locations of the street trees

  • Use the rasterstats package to calculate the NDVI values at the locations of the street trees.
  • Since these are point geometries, you can calculate either the median or the mean statistic (only one pixel will contain each point).
In [ ]:
#calculate NDVI values with zonal_stats
from rasterstats import zonal_stats
tree_stats = zonal_stats(tree, NDVI_city, affine=landsat.transform, stats=[ 'median'])
tree_stats

2.2.3 Plotting the results

Make two plots of the results:

  1. A histogram of the NDVI values, using matplotlib's hist function. Include a vertical line that marks the NDVI = 0 threshold
  2. A plot of the street tree points, colored by the NDVI value, using geopandas' plot function. Include the city limits boundary on your plot.

The figures should be clear and well-styled, with for example, labels for axes, legends, and clear color choices.

In [45]:
# Store the median value in the tree data frame
tree['median_NDVI'] = [s['median'] for s in tree_stats]
tree.head()
Out[45]:
objectid fcode geometry median_NDVI
0 1 3000 POINT (499541.269 4434698.265) 0.235337
1 2 3000 POINT (488932.471 4424093.158) 0.261535
2 3 3000 POINT (489039.214 4425985.827) 0.096769
3 4 3000 POINT (488993.171 4426088.005) 0.076630
4 5 3000 POINT (488943.113 4424599.478) 0.267952
In [46]:
#plot hist
fig, ax = plt.subplots(figsize=(10, 6))

ax.hist(tree["median_NDVI"],30, density=True, facecolor='g', alpha=0.75)

# Format the axes
ax.set_xlabel("Median NDVI")
ax.set_ylabel("Number of trees")
ax.set_axisbelow(True)
ax.grid(True, axis="y",alpha=0.5,linestyle="dotted")
ax.set_yticklabels([f"{yval:,.0f}" for yval in ax.get_yticks()] )
ax.set_title("Median NDVI of Trees in Philadelphia")

ax.spines['top'].set_visible(False)
ax.spines['right'].set_visible(False)

# Add the NDVI = 0 threshold
ax.axvline(x=0, c='k', linewidth=2)
ax.text(0, 3, " NDVI = 0", ha='left', fontsize=14);
F:\anaconda\envs\musa-550-fall-2020\lib\site-packages\ipykernel_launcher.py:11: UserWarning: FixedFormatter should only be used together with FixedLocator
  # This is added back by InteractiveShellApp.init_path()
In [47]:
import contextily as ctx
In [48]:
#A plot of the street tree points, colored by the NDVI values
# create the axes
fig, ax = plt.subplots(figsize=(14, 14))

# add background tracts
evict_3857.to_crs(tree.crs).plot(ax=ax, facecolor="none", edgecolor="white",alpha=0.75, linewidth=0.25)

# plot trees
tree.plot(ax=ax,column='median_NDVI', marker='.',legend=True)

# add the city limits
city_limits.to_crs(tree.crs).plot(ax=ax, edgecolor='white', linewidth=3, facecolor='none')

# plot the basemap underneath
ctx.add_basemap(ax=ax, crs=tree.crs, source=ctx.providers.CartoDB.DarkMatter)

# remove axis lines
ax.set_axis_off()
ax.set_title("Median NDVI of Trees in Philadelphia")
Out[48]:
Text(0.5, 1.0, 'Median NDVI of Trees in Philadelphia')